#include <map>
#include <mpi.h>
#include <omp.h>
#include <cmath>
#include <time.h>
#include <vector>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <iostream>
#include <string.h>
#include <windows.h>
#include <algorithm>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
using namespace std;

#define NUM_THREADS 8

const int maxsize = 3000;
const int maxrow = 40000;
const int numBasis = 40000;
int num;
map<int, int*>ans;
vector<int> tmpAns;
long long head, tail, freq;

fstream RowFile("被消元行.txt", ios::in | ios::out);
fstream BasisFile("消元子.txt", ios::in | ios::out);

int gRows[maxrow][maxsize];
int gBasis[numBasis][maxsize];
int answers[maxrow][maxsize];
map<int, int>firstToRow;

int ifBasis[numBasis] = { 0 };
int ifDone[maxrow] = { 0 };

void reset()
{
	memset(gRows, 0, sizeof(gRows));
	memset(gBasis, 0, sizeof(gBasis));
	memset(ifBasis, 0, sizeof(ifBasis));
	RowFile.close();
	BasisFile.close();
	RowFile.open("被消元行.txt", ios::in | ios::out);
	BasisFile.open("消元子.txt", ios::in | ios::out);
	ans.clear();
}

int readBasis()
{
	for (int i = 0; i < numBasis; i++) {
		if (BasisFile.eof()) {
			cout << "读取消元子" << i - 1 << "行" << endl;
			return i - 1;
		}
		string tmp;
		bool flag = false;
		int row = 0;
		getline(BasisFile, tmp);
		stringstream s(tmp);
		int pos;
		while (s >> pos) {
			if (!flag) {
				row = pos;
				flag = true;
				ifBasis[row] = 1;
			}
			int index = pos / 32;
			int offset = pos % 32;
			gBasis[row][index] = gBasis[row][index] | (1 << offset);
		}
		flag = false;
		row = 0;
	}
}

int readRowsFrom(int pos)
{
	if (RowFile.is_open())
		RowFile.close();
	RowFile.open("被消元行.txt", ios::in | ios::out);
	memset(gRows, 0, sizeof(gRows));   //重置为0
	string line;
	for (int i = 0; i < pos; i++) {       //读取pos前的无关行
		getline(RowFile, line);
	}
	for (int i = pos; i < pos + maxrow; i++) {
		int tmp;
		getline(RowFile, line);
		if (line.empty()) {
			cout << "读取被消元行 " << i << " 行" << endl;
			return i;   //返回读取的行数
		}
		bool flag = false;
		stringstream s(line);
		while (s >> tmp) {
			int index = tmp / 32;
			int offset = tmp % 32;
			gRows[i - pos][index] = gRows[i - pos][index] | (1 << offset);
			flag = true;
		}
	}
	cout << "read max rows" << endl;
	return -1;  //成功读取maxrow行
}

int findfirst(int row)
{
	int first;
	for (int i = maxsize - 1; i >= 0; i--) {
		if (gRows[row][i] == 0)
			continue;
		else {
			int pos = i * 32;
			int offset = 0;
			for (int k = 31; k >= 0; k--) {
				if (gRows[row][i] & (1 << k)) {
					offset = k;
					break;
				}
			}
			first = pos + offset;
			return first;
		}
	}
	return -1;
}

int _findfirst(int row)
{
	int first;
	for (int i = maxsize - 1; i >= 0; i--) {
		if (answers[row][i] == 0)
			continue;
		else {
			int pos = i * 32;
			int offset = 0;
			for (int k = 31; k >= 0; k--) {
				if (answers[row][i] & (1 << k)) {
					offset = k;
					break;
				}
			}
			first = pos + offset;
			return first;
		}
	}
	return -1;
}

void writeResult(ofstream& out)
{
	for (auto it = ans.rbegin(); it != ans.rend(); it++) {
		int* result = it->second;
		int max = it->first / 32 + 1;
		for (int i = max; i >= 0; i--) {
			if (result[i] == 0)
				continue;
			int pos = i * 32;
			for (int k = 31; k >= 0; k--)
				if (result[i] & (1 << k))
					out << k + pos << " ";
		}
		out << endl;
	}
}

void grobner()
{
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);
	int num = (flag == -1) ? maxrow : flag;
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	for (int i = 0; i < num; i++) {
		while (findfirst(i) != -1) {     //存在首项
			int first = findfirst(i);      //first是首项
			if (ifBasis[first] == 1) {
				for (int j = 0; j < maxsize; j++)
					gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
			}
			else {
				for (int j = 0; j < maxsize; j++)
					gBasis[first][j] = gRows[i][j];
				ifBasis[first] = 1;
				ans.insert(pair<int, int*>(first, gBasis[first]));
				break;
			}
		}
	}
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "grobner" << (tail - head) * 1000 / freq << "ms" << endl;
}

void g_omp()
{
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);     //读取被消元行
	int t_id = omp_get_thread_num();
	int num = (flag == -1) ? maxrow : flag;
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
#pragma omp parallel num_threads(NUM_THREADS)
	{
#pragma omp for schedule(guided)
		for (int i = 0; i < num; i++) {
			while (findfirst(i) != -1) {     //存在首项
				int first = findfirst(i);      //first是首项
				if (ifBasis[first] == 1) {
					for (int j = 0; j < maxsize; j++) {
						gRows[i][j] = gRows[i][j] ^ gBasis[first][j];     //进行异或消元
					}
				}
				else {
#pragma omp critical
					if (ifBasis[first] == 0) {
						for (int j = 0; j < maxsize; j++) {
							gBasis[first][j] = gRows[i][j];
						}
						ifBasis[first] = 1;
						ans.insert(pair<int, int*>(first, gBasis[first]));
					}
				}
			}
		}
	}
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "g_omp" << (tail - head) * 1000 / freq << "ms" << endl;
}

void g_sse()
{
    int begin = 0;
    int flag;
    flag = readRowsFrom(begin);     //读取被消元行
    int num = (flag == -1) ? maxrow : flag;
    LARGE_INTEGER head, tail, freq;
    QueryPerformanceCounter(&head);
    for (int i = 0; i < num; i++) {
        while (findfirst(i) != -1) {
            int first = findfirst(i);
            if (ifBasis[first] == 1) {  //存在该消元子
                int j = 0;
                for (; j + 4 < maxsize; j += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*) &gRows[i][j]);
                    __m128i vj = _mm_loadu_si128((__m128i*) &gBasis[first][j]);
                    __m128i vx = _mm_xor_si128(vij, vj);
                    _mm_storeu_si128((__m128i*) &gRows[i][j], vx);
                }
                for (; j < maxsize; j++)
                    gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
            }
            else {
                int j = 0;
                for (; j + 4 < maxsize; j += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*) &gRows[i][j]);
                    _mm_storeu_si128((__m128i*) &gBasis[first][j], vij);
                }
                for (; j < maxsize; j++)
                    gBasis[first][j] = gRows[i][j];
                ifBasis[first] = 1;
                ans.insert(pair<int, int*>(first, gBasis[first]));
                break;
            }
        }
    }
    QueryPerformanceCounter(&tail);
    QueryPerformanceFrequency(&freq);
    cout << "g_sse" << (tail.QuadPart - head.QuadPart) * 1000.0 / freq.QuadPart << "ms" << endl;
}

void writeResult_mpi(ofstream& out)
{
	for (int j = 0; j < num; j++) {
		for (int i = maxsize - 1; i >= 0; i--) {
			if (answers[j][i] == 0)
				continue;
			int pos = i * 32;
			//int offset = 0;
			for (int k = 31; k >= 0; k--) {
				if (answers[j][i] & (1 << k)) {
					out << k + pos << " ";
				}
			}
		}
		out << endl;
	}
}

void g_mpi(int argc, char* argv[])
{
	int flag;
	double start = 0, end_time = 0;
	MPI_Init(&argc, &argv);
	int total = 0, rank = 0, i = 0, j = 0, begin = 0, end = 0;
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &total);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	if (rank == 0) {
		flag = readRowsFrom(0);     //读取被消元行
		num = (flag == -1) ? maxrow : flag;
		begin = rank * num / total;
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = 1; i < total; i++) {
			MPI_Send(&num, 1, MPI_INT, i, 0, MPI_COMM_WORLD);//0是被消元行行数
			int b = i * (num / total);
			int e = (i == total - 1) ? num : (i + 1) * (num / total);
			for (j = b; j < e; j++) {
				MPI_Send(&gRows[j][0], maxsize, MPI_INT, i, 1, MPI_COMM_WORLD);//1时被消元行数据
			}
		}
	}
	else {
		MPI_Recv(&num, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		begin = rank * (num / total);
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = begin; i < end; i++) {
			MPI_Recv(&gRows[i][0], maxsize, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
	start = MPI_Wtime();
	for (i = begin; i < end; i++) {
		int first = findfirst(i);
		while (first != -1) {
			if (ifBasis[first] == 1) {  //存在首项为first消元子
				for (j = 0; j < maxsize; j++)
					gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
				first = findfirst(i);
			}
			else {   //升级为消元子
				tmpAns.push_back(first);
				if (rank == 0) {
					for (j = 0; j < maxsize; j++) {
						gBasis[first][j] = gRows[i][j];
						answers[i][j] = gRows[i][j];
					}
					ifBasis[first] = 1;  //仅仅将0号进程消元到底
				} break;
			}
		}
		if (first == -1)
			tmpAns.push_back(-1);
	}
	for (i = 0; i < rank; i++) {
		int b = i * (num / total);
		int e = b + num / total;
		for (j = b; j < e; j++) {
			MPI_Recv(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD, &status);//接收来自进程i的消元结果，可能作为之后的消元子
			int first = _findfirst(j);
			firstToRow.insert(pair<int, int>(first, j));//记录下首项信息
		}
		for (j = begin; j < end; j++) { 
			int first = tmpAns.at(j - begin);
			if (first == -1)
				continue;
			while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
				if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
					int row = firstToRow.find(first)->second;
					for (int k = 0; k < maxsize; k++)
						gRows[j][k] = gRows[j][k] ^ answers[row][k];
					first = findfirst(i);
				}
				if (first == -1)
					break;
				if (ifBasis[first] == 1) {
					for (int k = 0; k < maxsize; k++)
						gRows[j][k] = gRows[j][k] ^ gBasis[first][k];
					first = findfirst(i);
				}
			}
		}
	}
	if (rank != 0) {
		for (i = begin; i < end; i++) {
			int first = findfirst(i);
			if (first == -1)
				continue;
			while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
				if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
					int row = firstToRow.find(first)->second;
					for (int k = 0; k < maxsize; k++)
						gRows[i][k] = gRows[i][k] ^ answers[row][k];
					first = findfirst(i);
				}

				if (first == -1)
					break;
				if (ifBasis[first] == 1) {
					for (int k = 0; k < maxsize; k++)
						gRows[i][k] = gRows[i][k] ^ gBasis[first][k];
					first = findfirst(i);
				}
			}
			for (j = 0; j < maxsize; j++) {
				gBasis[first][j] = gRows[i][j];
				answers[i][j] = gRows[i][j];  //自身进程的消元结果不会加入firstToRow
			}
			ifBasis[first] = 1;
		}
	}
	for (i = rank + 1; i < total; i++) {
		for (j = begin; j < end; j++)
			MPI_Send(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD);
	}

	if (rank == total - 1) {
		end_time = MPI_Wtime();
		cout << "g_mpi" << 1000 * (end_time - start) << "ms" << endl;
		writeResult_MPI(out_mpi);
		out_mpi.close();
	}
	MPI_Finalize();
}

void g_mpi_omp(int argc, char* argv[])
{
	int flag;
	double start = 0;
	double end_time = 0;
	MPI_Init(&argc, &argv);
	int total = 0, rank = 0, i = 0, j = 0, begin = 0, end = 0;
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &total);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	if (rank == 0) {
		flag = readRowsFrom(0);     //读取被消元行
		num = (flag == -1) ? maxrow : flag;
		begin = rank * num / total;
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = 1; i < total; i++) {
			MPI_Send(&num, 1, MPI_INT, i, 0, MPI_COMM_WORLD);//0是被消元行行数
			int b = i * (num / total);
			int e = (i == total - 1) ? num : (i + 1) * (num / total);
			for (j = b; j < e; j++) {
				MPI_Send(&gRows[j][0], maxsize, MPI_INT, i, 1, MPI_COMM_WORLD);//1时被消元行数据
			}
		}

	}
	else {
		MPI_Recv(&num, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		begin = rank * (num / total);
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = begin; i < end; i++) {
			MPI_Recv(&gRows[i][0], maxsize, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);
	start = MPI_Wtime();
#pragma omp parallel  num_threads(NUM_THREADS),private(i,j)
#pragma omp for ordered schedule(guided)
	for (i = begin; i < end; i++) {
		int first = findfirst(i);

		while (first != -1) {
			if (ifBasis[first] == 1) {  //存在首项为first消元子
				for (j = 0; j < maxsize; j++) {
					gRows[i][j] = gRows[i][j] ^ gBasis[first][j];     //进行异或消元
				}
				first = findfirst(i);
			}
			else {
				if (rank == 0) {
#pragma omp ordered
					for (j = 0; j < maxsize; j++) {
						gBasis[first][j] = gRows[i][j];
						answers[i][j] = gRows[i][j];
					}
					ifBasis[first] = 1;  //仅仅将0号进程消元到底
				}
				first = -1;
			}
		}
	}
	for (i = 0; i < rank; i++) {
		int b = i * (num / total);
		int e = b + num / total;
		for (j = b; j < e; j++) {
			MPI_Recv(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD, &status);//接收来自进程i的消元结果，可能作为之后的消元子
			int first = _findfirst(j);
			firstToRow.insert(pair<int, int>(first, j));//记录下首项信息
		}
#pragma omp for schedule(guided)
		for (j = begin; j < end; j++) {
			int first = findfirst(j);
			while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
				if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
					int row = firstToRow.find(first)->second;
					for (int k = 0; k < maxsize; k++) 
						gRows[j][k] = gRows[j][k] ^ answers[row][k];
					first = findfirst(i);
				}

				if (ifBasis[first] == 1 && first != -1) {
					for (int k = 0; k < maxsize; k++) 
						gRows[j][k] = gRows[j][k] ^ gBasis[first][k];
					first = findfirst(i);
				}
			}
		}
	}
	if (rank != 0) {
		for (i = begin; i < end; i++) {
			int first = findfirst(i);
			while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
				if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
					int row = firstToRow.find(first)->second;
					for (int k = 0; k < maxsize; k++)
						gRows[i][k] = gRows[i][k] ^ answers[row][k];
					first = findfirst(i);
				}
				if (ifBasis[first] == 1 && first != -1) {
					for (int k = 0; k < maxsize; k++)
						gRows[i][k] = gRows[i][k] ^ gBasis[first][k];
					first = findfirst(i);
				}
			}
			for (j = 0; j < maxsize; j++) {
				gBasis[first][j] = gRows[i][j];
				answers[i][j] = gRows[i][j];  //自身进程的消元结果不会加入firstToRow
			}
			ifBasis[first] = 1;
		}
	}
	for (i = rank + 1; i < total; i++) {
		for (j = begin; j < end; j++) {
			MPI_Send(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD);
		}
	}
	if (rank == total - 1) {
		end_time = MPI_Wtime();
		cout << "g_mpi_omp " << 1000 * (end_time - start) << "ms" << endl;
		writeResult_MPI(out_mpi);
		out_mpi.close();
	}
	MPI_Finalize();
}

void g_mpi_sse(int argc, char* argv[])
{
	int flag;
	double start = 0, end_time = 0;
	MPI_Init(&argc, &argv);
	int total = 0, rank = 0, i = 0, j = 0, begin = 0, end = 0;
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &total);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	if (rank == 0) {
		flag = readRowsFrom(0);     //读取被消元行
		num = (flag == -1) ? maxrow : flag;
		begin = rank * num / total;
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = 1; i < total; i++) {
			MPI_Send(&num, 1, MPI_INT, i, 0, MPI_COMM_WORLD);//0是被消元行行数
			int b = i * (num / total);
			int e = (i == total - 1) ? num : (i + 1) * (num / total);
			for (j = b; j < e; j++) {
				MPI_Send(&gRows[j][0], maxsize, MPI_INT, i, 1, MPI_COMM_WORLD);//1时被消元行数据
			}
		}
	}
	else {
		MPI_Recv(&num, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		begin = rank * (num / total);
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = begin; i < end; i++) {
			MPI_Recv(&gRows[i][0], maxsize, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);
	start = MPI_Wtime();

	for (i = begin; i < end; i++) {
        int first = findfirst(i);
        while (first != -1) { // 未消元完毕，存在首项
            if (ifBasis[first] == 1) { // 存在首项为first消元子
                for (j = 0; j + 4 < maxsize; j += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*)&gRows[i][j]);
                    __m128i vj = _mm_loadu_si128((__m128i*)&gBasis[first][j]);
                    __m128i vx = _mm_xor_si128(vij, vj);
                    _mm_storeu_si128((__m128i*)&gRows[i][j], vx);
                }
                for (; j < maxsize; j++) {
                    gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
                }
                first = findfirst(i);
            } else {
                if (rank == 0) {
                    for (j = 0; j + 4 < maxsize; j += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*)&gRows[i][j]);
                        _mm_storeu_si128((__m128i*)&gBasis[first][j], vij);
                        _mm_storeu_si128((__m128i*)&answers[i][j], vij);
                    }
                    for (; j < maxsize; j++) {
                        gBasis[first][j] = gRows[i][j];
                        answers[i][j] = gRows[i][j];
                    }
                    ifBasis[first] = 1;
                }
                break;
            }
        }
    }


    for (i = 0; i < rank; i++) {
        int b = i * (num / total);
        int e = b + num / total;
        for (j = b; j < e; j++) {
            MPI_Recv(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD, &status);
            int first = _findfirst(j);
            firstToRow.insert(pair<int, int>(first, j));
        }
        for (j = begin; j < end; j++) {
            int first = findfirst(j);
            while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
                if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
                    int row = firstToRow.find(first)->second;
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & answers[row][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[j][k] = gRows[j][k] ^ answers[row][k];
                    }
                    first = findfirst(i);
                }
                if (ifBasis[first] == 1) {
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[j][k] = gRows[j][k] ^ gBasis[first][k];
                    }
                    first = findfirst(i);
                }
            }
        }
    }

    if (rank != 0) {
        for (i = begin; i < end; i++) {
            int first = findfirst(i);
            if (first == -1)
                continue;
            while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
                if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
                    int row = firstToRow.find(first)->second;
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & answers[row][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[i][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[i][k] = gRows[i][k] ^ answers[row][k];
                    }
                    first = findfirst(i);
                }
                if (first == -1)
                    break;
                if (ifBasis[first] == 1) {
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[i][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[i][k] = gRows[i][k] ^ gBasis[first][k];
                    }
                    first = findfirst(i);
                }
            }
            if (first != -1) {
                for (j = 0; j < maxsize; j++) {
                    gBasis[first][j] = gRows[i][j];
                    answers[i][j] = gRows[i][j];  //自身进程的消元结果不会加入firstToRow
                }
                ifBasis[first] = 1;
            }
        }
    }

	for (i = rank + 1; i < total; i++)
		for (j = begin; j < end; j++)
			MPI_Send(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD);

	if (rank == total - 1) {
		end_time = MPI_Wtime();
		cout << "g_mpi_avx " << 1000 * (end_time - start) << "ms" << endl;
		writeResult_MPI(out_mpi);
		out_mpi.close();
	}
	MPI_Finalize();
}

void g_mpi_sse_omp(int argc, char* argv[])
{
	int flag;
	double start = 0;
	double end_time = 0;
	MPI_Init(&argc, &argv);
	int total = 0, rank = 0, i = 0, j = 0, begin = 0, end = 0;
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &total);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	if (rank == 0) {
		flag = readRowsFrom(0);     //读取被消元行
		num = (flag == -1) ? maxrow : flag;
		begin = rank * num / total;
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = 1; i < total; i++) {
			MPI_Send(&num, 1, MPI_INT, i, 0, MPI_COMM_WORLD);//0是被消元行行数
			int b = i * (num / total);
			int e = (i == total - 1) ? num : (i + 1) * (num / total);
			for (j = b; j < e; j++) {
				MPI_Send(&gRows[j][0], maxsize, MPI_INT, i, 1, MPI_COMM_WORLD);//1时被消元行数据
			}
		}

	}
	else {
		MPI_Recv(&num, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
		begin = rank * (num / total);
		end = (rank == total - 1) ? num : (rank + 1) * (num / total);
		for (i = begin; i < end; i++) {
			MPI_Recv(&gRows[i][0], maxsize, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
		}
	}


	MPI_Barrier(MPI_COMM_WORLD);  //此时每个进程都拿到了数据
	start = MPI_Wtime();
#pragma omp parallel  num_threads(NUM_THREADS),private(i,j)
#pragma omp for ordered schedule(guided)
    for (i = begin; i < end; i++) {
        int first = findfirst(i);
        while (first != -1) {
            if (ifBasis[first] == 1) {  //存在首项为first消元子
                for (j = 0; j + 4 < maxsize; j += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][j]);
                    __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][j]);
                    __m128i vx = _mm_xor_si128(vij, vj);
                    _mm_storeu_si128((__m128i*) & gRows[i][j], vx);
                }
                for (; j < maxsize; j++) {
                    gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
                }
                first = findfirst(i);
            }
            else {
    #pragma omp ordered
                if (rank == 0) {
                    while (ifBasis[first] == 1) {
                        for (j = 0; j + 4 < maxsize; j += 4) {
                            __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][j]);
                            __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][j]);
                            __m128i vx = _mm_xor_si128(vij, vj);
                            _mm_storeu_si128((__m128i*) & gRows[i][j], vx);
                        }
                        for (; j < maxsize; j++)
                            gRows[i][j] = gRows[i][j] ^ gBasis[first][j];
                        first = findfirst(i);
                    }
                    if (first != -1) {
                        for (j = 0; j + 4 < maxsize; j += 4) {
                            __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][j]);
                            _mm_storeu_si128((__m128i*) & gBasis[first][j], vij);
                            _mm_storeu_si128((__m128i*) & answers[i][j], vij);
                        }
                        for (; j < maxsize; j++) {
                            gBasis[first][j] = gRows[i][j];
                            answers[i][j] = gRows[i][j];
                        }
                        ifBasis[first] = 1;
                    }
                }
                first = -1;
            }
        }
    }
  
	for (i = 0; i < rank; i++) {
		int b = i * (num / total);
		int e = b + num / total;
		for (j = b; j < e; j++) {
			MPI_Recv(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD, &status);//接收来自进程i的消元结果，可能作为之后的消元子
			int first = _findfirst(j);
			firstToRow.insert(pair<int, int>(first, j));//记录下首项信息
		}
        #pragma omp for schedule(guided)
        for (j = begin; j < end; j++) {
            int first = findfirst(j);
            while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
                if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
                    int row = firstToRow.find(first)->second;
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & answers[row][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[j][k] = gRows[j][k] ^ answers[row][k];
                    }
                    first = findfirst(i);
                }
                if (ifBasis[first] == 1) {
                    int k = 0;
                    for (k = 0; k + 4 < maxsize; k += 4) {
                        __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                        __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][k]);
                        __m128i vx = _mm_xor_si128(vij, vj);
                        _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                    }
                    for (; k < maxsize; k++) {
                        gRows[j][k] = gRows[j][k] ^ gBasis[first][k];
                    }
                    first = findfirst(i);
                }
            }
        }
        int first = findfirst(j);
        while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
            if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
                int row = firstToRow.find(first)->second;
                int k = 0;
                for (k = 0; k + 4 < maxsize; k += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                    __m128i vj = _mm_loadu_si128((__m128i*) & answers[row][k]);
                    __m128i vx = _mm_xor_si128(vij, vj);
                    _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                }
                for (; k < maxsize; k++) {
                    gRows[j][k] = gRows[j][k] ^ answers[row][k];
                }
                first = findfirst(i);
            }
            if (ifBasis[first] == 1) {
                int k = 0;
                for (k = 0; k + 4 < maxsize; k += 4) {
                    __m128i vij = _mm_loadu_si128((__m128i*) & gRows[j][k]);
                    __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][k]);
                    __m128i vx = _mm_xor_si128(vij, vj);
                    _mm_storeu_si128((__m128i*) & gRows[j][k], vx);
                }
                for (; k < maxsize; k++) {
                    gRows[j][k] = gRows[j][k] ^ gBasis[first][k];
                }
                first = findfirst(i);
            }
        }
	}

    if (rank != 0) {
        for (i = begin; i < end; i++) {
            int first = findfirst(i);
            if (first != -1) {
                while ((firstToRow.find(first) != firstToRow.end() || ifBasis[first] == 1) && first != -1) {  //存在可消元项
                    if (firstToRow.find(first) != firstToRow.end()) {  //消元结果有消元子
                        int row = firstToRow.find(first)->second;
                        int k = 0;
                        for (k = 0; k + 4 < maxsize; k += 4) {
                            __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][k]);
                            __m128i vj = _mm_loadu_si128((__m128i*) & answers[row][k]);
                            __m128i vx = _mm_xor_si128(vij, vj);
                            _mm_storeu_si128((__m128i*) & gRows[i][k], vx);
                        }
                        for (; k < maxsize; k++)
                            gRows[i][k] = gRows[i][k] ^ answers[row][k];
                        first = findfirst(i);
                    }
                    if (ifBasis[first] == 1) {
                        int k = 0;
                        for (k = 0; k + 4 < maxsize; k += 4) {
                            __m128i vij = _mm_loadu_si128((__m128i*) & gRows[i][k]);
                            __m128i vj = _mm_loadu_si128((__m128i*) & gBasis[first][k]);
                            __m128i vx = _mm_xor_si128(vij, vj);
                            _mm_storeu_si128((__m128i*) & gRows[i][k], vx);
                        }
                        for (; k < maxsize; k++)
                            gRows[i][k] = gRows[i][k] ^ gBasis[first][k];
                        first = findfirst(i);
                    }
                }
                if (first == -1) {
                    continue;
                }
                for (j = 0; j < maxsize; j++) {
                    gBasis[first][j] = gRows[i][j];
                    answers[i][j] = gRows[i][j];
                }
                ifBasis[first] = 1;
            }
        }
    }

	for (i = rank + 1; i < total; i++) {
		for (j = begin; j < end; j++)
			MPI_Send(&answers[j][0], maxsize, MPI_INT, i, 2, MPI_COMM_WORLD);
	}

	if (rank == total - 1) {
		end_time = MPI_Wtime();
		cout << "g_mpi_avx_omp " << 1000 * (end_time - start) << "ms" << endl;
		writeResult_MPI(out_mpi);
		out_mpi.close();
	}
	MPI_Finalize();
}

int main(int argc, char* argv[])
{
	QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
	readBasis();
	mpi(argc, argv);
}